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The  detailed  isothermal  electrochemical  model  for  a  lithium  ion  cell  comprises  of  10  coupled  partial 
differential  equations  (PDEs).  It  is  computationally  intensive  to  use  this  model  for  detailed  parameter 
estimation,  battery  packs,  battery  management  system  algorithms,  calendar  or  cycle  life  predictions  etc. 
This  work  describes  a  reduced  order  model  framework,  which  reduces  the  detailed  model  PDEs,  to 
a  manageable  set  of  ordinary  differential  equations  (ODEs),  which  can  be  used  for  the  above  applications. 
Volume  averaging  PDEs  yield  ODEs.  Hence  volume  averaging  is  the  basic  method  used  for  model 
reduction.  The  gradient  or  profile  information  lost  on  volume  averaging  is  recovered  through  profile 
based  approximations.  For  the  lowest  order  reduced  order  model  (ROM),  a  uniform  reaction  rate, 
quadratic  electrolyte  concentrations  and  quartic  solid  phase  concentrations  are  assumed.  In  this  ROM 
only  5  linear  ODEs  are  solved,  and  the  rest  are  nonlinear  algebraic  evaluations.  When  compared  with  the 
detailed  electrochemical  model  of  a  lithium  ion  cell,  this  ROM  gives  negligible  error  at  nominal  currents 
and  a  maximum  error  of  about  3%  at  high  currents,  for  a  typical  commercial  cell. 

©  2012  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Consider  a  lithium  ion  cell  discharging  current  through  a  load,  as 
in  Fig.  1.  Say  the  negative  electrode  is  graphite.  When  the  cell  is  in 
charged  state  lithium  exists  as  the  lithium  carbide  (LiCe)  in  the 
negative  electrode.  During  discharge  LiC6  breaks  down  releasing 
Li+  ions  into  the  electrolyte,  and  the  electrons  to  the  collector 
(through  the  solid  conducting  matrix  of  the  negative  electrode 
particles),  which  eventually  flow  into  the  external  circuit.  The 
electrons  travel  from  the  negative  electrode  to  the  positive  elec¬ 
trode  through  the  external  circuit,  while  the  Li+  ions  travel  in  the 
electrolyte  from  the  negative  electrode  to  the  positive  electrode. 
Say  the  positive  electrode  is  some  lithium  metal  oxide.  The  elec¬ 
trons  (flowing  through  the  external  circuit)  and  Li+  ions  (diffusing 
through  the  electrolyte)  combine  on  the  surface  of  the  positive 
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electrode  particles  to  form  neutral  Li  species,  which  then  interca¬ 
late  into  lithium  metal  oxide.  The  intercalation  reactions  occurring 
during  the  discharge  are  reversible  and  the  reverse  reactions  occur 
during  charge. 

Electrodes  contain  a  solid  phase  comprising  of  active  materials 
(and  other  additives),  and  a  liquid  electrolyte  phase.  Separator 
typically  has  some  solid  polymeric  membrane  and  a  liquid  elec¬ 
trolyte  phase.  Neglecting  the  solid  additives  and  inert  polymeric 
membranes,  there  are  totally  five  phases:  active  materials  in  the 
two  electrodes  and  electrolyte  phase  in  the  three  regions.  For  an 
isothermal  operation  of  a  single  cell,  drawing  the  mass  and  charge 
balances  for  these  five  phases,  the  detailed  electrochemical  model 
comprises  of  ten  coupled  partial  differential  equations  (PDEs),  Refs. 
[1,2],  It  is  computationally  intensive  (or  often  impossible)  to  use 
this  detailed  electrochemical  model  for 

•  Battery  packs,  since  they  have  hundreds  of  single  cells  in  series 
-  parallel  configurations. 

•  Battery  management  system  algorithms,  since  on-board 
computational  capabilities  are  nominal.  Currently  electrical 
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LiC6— >x  Li++Li1_xC6+x  e_  x  Li++  LiM02+x  e-  — >  Li1+xM02 

Fig.  1.  Discharge  of  a  lithium  ion  cell. 


circuit  based  models  are  used  for  battery  management 
systems,  due  to  their  minimal  mathematical  complexity  (first 
order  ordinary  differential  equations  based  mathematical 
structure). 

•  Calendar  or  cycle  life  predictions,  since  the  model  needs  to 
predict  cell  behavior  for  real  time  periods  as  long  as  a  few 
months. 

•  Detailed  parameters  estimation,  since  systematic  optimization 
of  parameters  embedded  within  ten  coupled  PDEs  is 
formidable. 

•  inclusion  of  complex  degradation  mechanisms,  since  they  add  to 
the  numerical  stiffness. 

The  aim  of  this  work  is  to  develop  a  reduced  order  model 
framework,  which  reduces  the  detailed  electrochemical  model 
comprising  of  ten  PDEs  for  the  isothermal  behavior  of  a  single  cell, 
to  a  manageable  set  of  ordinary  differential  equations  (ODEs), 
which  can  be  used  for  the  above  applications.  It  is  preferred  that  the 
reduced  order  model  framework  has  the  following  mathematical 
properties: 

•  It  has  minimal  computational  requirements  so  that  it  can  be 
implemented  in  on-board  battery  management  systems.  It  is 
numerically  stable  and  easily  optimized  during  parameter 
estimations.  These  properties  are  ensured  if  the  model  equa¬ 
tions  are  first  order  ODEs  (linear  or  nonlinear,  could  be 
coupled).  Volume  averaging  PDEs  yield  ODEs.  Hence  volume 
averaging  is  the  basic  methodology  used  in  this  work  for 
physics  based  model  reduction. 

•  The  model  predictions  should  be  grid  or  system  size  inde¬ 
pendent.  For  example,  one  could  discretize  the  governing  PDEs, 
with  say  10  or  20  nodes  for  each  electrode  or  separator  region, 
as  in  the  Method  of  lines.  Model  predictions  in  such  a  scenario 
will  often  be  grid  size  dependent,  and  will  give  higher  errors  if 
the  electrode  or  separator  regions  get  thicker.  Thus,  ODEs 
should  be  solved  not  for  the  discretized  original  PDEs,  but  for 
physically  relevant  internal  variables,  which  parameterize  the 
internal  field  profiles.  For  example,  in  profile  based  approxi¬ 
mations,  interfacial  concentrations  are  solved,  and  the 
concentration  profiles  in  the  individual  regions  are  expressed 
in  terms  of  these  interfacial  concentrations.  A  brief  literature 
review  of  profile  based  approximations  is  given  below. 

•  It  has  systematic  mathematical  approximations  rather  than 
phenomenological  approximations.  For  example,  a  single 
particle  model  Ref.  [3]  neglects  electrolyte  potential  and  elec¬ 
trolyte  concentration  variations.  Hence,  a  single  particle  model 


cannot  work  at  high  currents.  If  a  reduced  order  model  is, 
however,  based  on  systematic  mathematical  approximations, 
relaxing  those  approximations  one  can  derive  higher  order 
models,  which  work  at  all  scenarios.  For  example,  the  current 
work  approximates  the  local  reaction  rates  in  an  electrode  to  be 
the  volume  average  reaction  rates.  This  approximation  can  be 
relaxed  to  allow  spatial  variation  of  reaction  rates,  will  be  re¬ 
ported  in  future. 

Solution  of  partial  differential  equations  through  profile  based 
approximations  was  developed  by  Pohlhausen  (1921 )  for  a  problem 
in  laminar  boundary  layer  theory,  Ref.  [4,  p.  41  ].  He  used  a  quartic 
velocity  profile  to  satisfy  the  momentum  integral  equation  and  the 
relevant  boundary  conditions.  Glueckauf  (1955)  [5]  used  a  para¬ 
bolic  concentration  profile  to  describe  the  diffusion  in  a  sphere,  in 
the  context  of  chromatography.  Ref.  [6]  used  the  parabolic 
approximation  for  diffusion  in  a  sphere  to  describe  the  kinetics  of 
adsorption  in  a  fixed  bed  of  spheres.  Refs.  [7,8],  for  example, 
describe  higher  order  polynomial  approximations  for  diffusion  in 
a  sphere.  Ref.  [9]  extended  these  methods  to  reaction  and  diffusion 
in  spherical  and  slab  geometries.  Ref.  [10]  reports  parabolic  and 
quartic  profile  based  approximations  for  diffusion  in  active  material 
spheres  in  the  electrodes  of  a  lithium  ion  cell.  This  work  extends  the 
profile  based  approximation  for  the  rectangular  electrolyte 
concentration  fields  in  electrode  and  separator  regions  as  well. 

Sections  2—5  present  the  detailed  field  equations  and  derive  the 
reduced  order  equations  for  solid  phase  current  balances,  electrolyte 
phase  mass  balances,  total  current  balances  and  solid  phase  mass 
balances,  respectively.  These  equations  are  the  core  conservation 
equations  governing  the  isothermal  behavior  of  a  single  cell.  The 
results  from  these  sections  feed  into  the  Butler— Volmer  kinetics,  in 
Section  6,  for  the  electrochemical  reactions  happening  on  the 
surface  of  the  active  material  spheres  in  the  electrode  regions.  The 
reduced  order  model  is  derived  across  Sections  2—6.  Often  the  deri¬ 
vation  sequence  is  different  from  the  computational  sequence.  Hence, 
the  reduced  order  model  results  are  collated  in  an  algorithmic  order  in 
Section  7.  The  cell  voltage  predictions  by  the  reduced  order  model  are 
compared  with  those  from  the  detailed  electrochemical  model  for 
constant  current  and  pulse  protocols  in  Section  8. 

2.  Solid  phase  current  balance  equations 

2.1.  Field  equations 

Ohm’s  law  is  typically  written  for  a  one-dimensional  conductor 
as  /  =  V/R,  where  I  is  the  current  in  Amperes  (A),  V  the  voltage  in 
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Volts  (V)  and  R  the  resistance  in  Ohms  (Q).  It  is  generalized  to 
a  three-dimensional  continuum  as  I  =  E/p,  where  I  is  the  current 
density  in  A  m-2,  E  is  the  electric  field  in  V/m,  and  p  the  resistivity  in 
Qm.  Conductivity,  a,  is  the  reciprocal  of  resistivity,  with  units 
Siemens  m  ',  where  1  Siemens  =  1  Q_1.  Thus,  Ohm’s  law  for 
a  continuum  is  /  =  -<rV$.  Ohm’s  law  describes  the  relationship 
between  current  density  and  potential  field  in  a  conductor.  In  the 
cell  electrodes,  the  electrochemical  reactions  generate  or  consume 
electrons.  Hence,  the  gradient  of  current  density  within  the  elec¬ 
trodes  is  proportional  to  electrochemical  reaction  rate.  These 
equations  are  the  solid  phase  current  balance  equations  in  the  two 
electrodes,  given  below  along  with  their  boundary  conditions.  At 
the  two  collector  ends,  solid  current  density  equals  the  total  current 
density,  computed  as  external  current  per  unit  electrode  cross 
sectional  area.  At  the  separator  interfaces,  solid  current  is  zero,  so 
that  within  the  separator  the  current  is  fully  ionic. 

2.1.1.  Negative  electrode  region 


Boundary  conditions  :  -  <7ln— ! -  =  I,  and  (2) 

dx  lx=0 

d&i  I 

=  0.  (3) 


2.1.2.  Positive  electrode  region 


r,  I 


Applying  the  boundary  conditions  (2)  and  (3)  as  the  lower  and 
upper  limits  of  integration  in  equation  (10)  give 


Similarly  volume  averaging  equation  (4)  gives 

I  =  -aplpF(jp).  (12) 

Equating  equations  (11)  and  (12)  and  rearranging  gives 

anln(jn)(t)  +  aplp(jp)(t)  =  0.  (13) 

In  a  closed  circuit,  the  electrons  generated  in  an  electrode  are 
consumed  in  the  other  electrode,  with  no  loss  of  electrons  in  the 
external  circuit.  Equation  (13)  accounts  for  this  overall  electro¬ 
neutrality. 

Equations  (11)  and  (12)  can  be  rearranged  to  get  an  expression 
for  the  average  reactions  rates  in  the  electrodes  as 

(in)  (t)  =  ^  and  (14) 


*(-a'pisr)  = 


-aPFjp, 


Boundary  condition  :  -  <rlp— —  =  0,  and  (5) 

°x  lx=(„+fs 

Q0.  I 

.=L  W 


2.2.  Volume  averaged  field  equations 


Volume  average  of  a  function  fix,t)  in  the  electrode  regions  is 
defined  as  follows: 


(f)n(t) 


Ain 


(7) 


<'5> 

Equations  (14)  and  (15)  relate  the  volume  average  reaction  rates  in 
the  electrodes  to  the  external  current,  at  any  time.  These  results  are 
exact,  since  volume  averaging  is  an  exact  operation  and  not  an 
approximation.  Of  course,  a  volume  averaged  field  equation  has  lost 
information  about  the  profile  or  gradient  of  a  quantity.  The  profile 
or  gradient  information  will  be  recovered  from  the  field  equation 
by  profile  approximation,  as  shown  below. 

For  a  constant  current  discharge,  the  volume  average  reaction 
rates  in  the  electrodes,  equations  (14)  and  (15),  are  constant.  If  by 
convention,  discharge  current  has  a  positive  sign,  then  the  volume 
average  reaction  rate  in  the  negative  electrode  is  positive,  while 
that  in  the  positive  electrode  is  negative.  These  signs  are  consistent 
with  the  fact  that  during  discharge  Li+  ions  are  released  in  the 
negative  electrode  region  and  consumed  in  the  positive  electrode 
region. 


(f}P(t)=;^  J  fix,t)A  dx 


2.3.  Profile  approximations 

(8)  In  the  uniform  reaction  rate  approximation,  the  local  reaction  rate 

is  assumed  to  be  the  volume  average  reaction  rate,  at  any  time  i.e. 


For  example  the  volume  average  reaction  rate  in  the  negative  jn(x,  t)  ~  (jn)(t)  =  ,  and 

electrode  is  0,1  n 


(16) 


=  J  jn(x,t)dx.  (9) 

x  =  0 

A  field  equation  is  volume  averaged  by  integrating  it  over  its 
domain,  and  normalizing  with  the  domain  volume.  For  example, 
equation  (1)  is  integrated  as  follows: 


Jp(x,t)  =  (/p)(t)  =  (17) 

The  errors  in  these  approximations  vanish  at  the  volume  averaged 
level. 

Any  spatio-temporal  quantity  can  be  written  as  a  sum  of  its  time 
varying  volume  average  and  a  spatio-temporal  deviation  term,  for 
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example  the  negative  electrode  reaction  rate  can  always  be  equality  across  interfaces.  At  the  collector  ends  zero  flux 
written  as  boundary  conditions  apply. 


jn(X,t)  =  (j„)(t)+j*(X,t). 


(18)  3.1.1.  Negative  electrode  region 


In  the  leading  order  uniform  reaction  rate  approximation,  the 
deviation  term  is  neglected.  This  work  explores  the  limits  of 
applicability  of  the  leading  order  volume  average  reaction  rate 
term. 

With  the  uniform  reaction  rate  approximation,  the  negative 
electrode  solid  current  balance  (1)  becomes 


(19) 


with  boundary  conditions 


The  above  equation  is  integrated  with  respect  to  x  to  get 


C2(ln,t)  =  C2in, 


(20) 


~D2n 


0C2 

0X 


=  ton- 


whereat)  is  some  arbitrary  function.  Using  the  boundary  condition  3.1.2.  Separator  region 
(3)  gives 


0  =  la+f(t )■ 


Eliminating /(t)  in  equation  (20)  using  equation  (21 )  gives  the  solid  with  boundary  conditions 

current  in  negative  electrode,  under  uniform  reaction  rate 

approximation,  as  c2(ln,t)  =  c2in, 


(22) 


=  ton, 


Note  that  equation  (22)  satisfies  the  other  boundary  condition  (2). 
Similarly  the  solid  current  in  the  positive  electrode,  under  uniform 
reaction  approximation  is 


=^[X-(/n  +  /s)] 


(23) 


C2(ln  +  lS,t)  =  C2i  p, 


and  -D2s 


0C2 

sdx 


\x=k+ls 


=  top- 


Thus  uniform  reaction  rate  approximation  implies  a  linear  solid 
current  profile  in  the  electrodes  i.e.  the  solid  current  density  line¬ 
arly  varies  from  the  total  current  density  at  the  collector  end  to  zero 
at  the  separator  end.  Hence,  a  test  of  the  solid  current  linearity  is 
a  test  of  uniform  reaction  rate  approximation;  the  numerical  results 
are  presented  in  supplementary  material  SI. 


3.1.3.  Positive  electrode  region 

£2p  if  =  i  (D2pSf) +  Qp(1  ~ t+)jp’ 

with  boundary  conditions 


3.  Electrolyte  phase  mass  balance  equations 


(24) 

(25) 

(26) 

(27) 

(28) 

(29) 

(30) 

(31) 

(32) 


(33) 

(34) 


Ref.  [10]  reports  profile  based  approximations  for  diffusion  in 
spherical  active  material  particles  in  the  electrodes.  This  work 
extends  that  approach  for  the  rectangular  electrolyte  concentration 
fields  in  the  electrode  and  separator  regions. 

3.1.  Field  equations 


-D2 


9C2 
3  dx 


lx=!„+/s 


=  top. 


and  -D2p^  =  °- 


(35) 

(36) 


The  electrolyte  diffusion  equations  form  the  electrolyte  phase 
mass  balances.  The  reaction  rate  terms  (described  in  the  previous 
section)  form  source  or  sink  terms  in  the  electrode  regions,  while 
electrolyte  accumulates  and  diffuses  through  the  separator.  Elec¬ 
trolyte  concentration  field  is  continuous  through  the  three  regions: 
negative  electrode,  separator  and  positive  electrode.  Hence,  the 
concentration  and  flux  continuities  hold  at  the  electrode-separator 
interfaces.  The  interfacial  concentrations  (c2m.C2ip)  and  fluxes 
(ton.top)  are  identified  as  separate  variables,  so  that  it  is  easy  to 
impose  them  in  the  individual  regions,  and  maintain  their 


Across  the  three  regions,  the  initial  condition  is 

c2(x,0)  =  c20.  (37) 

3.2.  Volume  averaged  field  equations 

The  volume  average  definitions  for  a  function  f{x,t)  in  the 
negative  and  positive  electrode  regions  are  as  defined  earlier  in 
equations  (7)  and  (8).  The  corresponding  definition  in  the  separator 
region  is 
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ln+ls  /n+/s 

(f)s(t)  =  ^  J  f(X,t)A  dx  J  /(x,t)dx.  (38) 

The  electrolyte  diffusion  equation  in  the  negative  electrode 
region,  equation  (24)  is  volume  averaged  as  follows: 


r,  J  h-ir  -  s  + -  «-]*-  <39> 


In  the  left  hand  side  of  equation  (39),  the  volume  average  of  the 
time  derivative  becomes  the  time  derivative  of  the  volume  average 


ms-¥  /  T&  =  «4  (s  /  C2dx)  -  '“TT-  (40> 

The  diffusive  flux  gradient  term  in  the  right  hand  side  of  equation 
(39)  gets  volume  averaged,  using  the  boundary  conditions  (25)  and 
(27)  as  follows: 


acA'" 

0*  /  x=0 


(41) 


The  volume  average  of  the  reaction  rate  term  in  the  right  hand  side 
of  equation  (39)  is  straight  forward: 


(an(l  -  t+)jn)n  =  /n  J  a„(l  -  t+)jn  dx  =  an(l  -  t+)(jn). 

x=0 

(42) 

Putting  together  the  results  from  equations  (40)— (42),  and  multi¬ 
plying  by  the  electrode  thickness,  the  volume  averaged  electrolyte 
diffusion  equation  in  the  negative  electrode  region  is 

inC2n%^  =  -q2,n  +  (1  -  t+)anln{jB).  (43) 

In  the  separator  region,  the  diffusive  flux  gradient  term  in  the 
right  hand  side  of  equation  (28)  gets  volume  averaged,  using  the 
boundary  conditions  (30)  and  (32)  as  follows: 


G«)>. 


_  ton  top 

~  %  '  is  ’ 

(44) 


go* 


H(d- 


0C2V 

3  0* /*=!„,  I, 


top 

lp 

(46) 


Volume  averaging  of  the  accumulation  and  reaction  rate  terms  in 
equation  (33)  is  as  before.  Finally,  the  volume  averaged  electrolyte 
diffusion  equation  in  the  positive  electrode  region  is 

d(C2)n  /.  \ 

h£2p — —  =  top  +  (1  -  1 1  )ap/p{jpy.  (47) 

The  evolution  equations  for  the  average  electrolyte  concentra¬ 
tions  in  the  three  regions,  equations  (43),  (45)  and  (47),  are  not 
independent  of  each  other.  During  discharge,  in  the  negative  elec¬ 
trode,  when  an  Li+  ion  gets  into  the  electrolyte  an  electron  gets 
released  into  the  outer  circuit.  This  electron  recombines  with  some 
other  Li+  in  the  electrolyte  within  the  positive  electrode  region, 
almost  instantaneously,  neglecting  electronic  transients.  Thus,  if 
a  Li+  ion  is  released  in  the  negative  electrode  region,  some  other  Li+ 
ion  gets  consumed  in  the  positive  electrode  region.  Hence,  an 
overall  electrolyte  balance  holds  at  any  time 

fnr2n(C2)n+isC2s(c2)s+fp£2p(C2)p  =  (inr2n+kr2s+fpE2p)c20!  (48) 


where  c2 o  is  the  initial  electrolyte  concentration,  as  in  the  initial 
condition  (37).  Differentiating  this  equation  with  respect  to  time 
gives 


(49) 


Adding  up  equations  (43),  (45)  and  (47),  and  using  the  electro¬ 
neutrality  condition  (13),  it  is  easily  seen  that  (49)  is  satisfied. 

In  the  volume  averaged  electrolyte  diffusion  equations  (43),  (45) 
and  (47),  the  interfacial  fluxes  appear  explicitly,  hence  they  are 
identified  as  physically  relevant  internal  variables  characterizing 
the  cell  behavior.  As  noted  earlier,  these  volume  averaged  results 
are  exact  and  not  approximations,  even  though  they  entail  loss  of 
concentration  profile  or  gradient  information. 


3.3.  Profile  approximations 
3.3. 1.  Negative  electrode  region 

While  deriving  the  volume  averaged  equation  for  the  negative 
electrode  region,  the  volume  average  of  the  diffusive  flux  gradient 
term  in  the  electrolyte  diffusion  equation  is  given  by  the  exact 
result,  equation  (41 ).  To  approximate  the  concentration  profile  in 
the  negative  electrode  region,  it  is  assumed  that 


Volume  averaging  of  the  accumulation  term  in  equation  (28)  is  as 
before.  Finally,  the  volume  averaged  electrolyte  diffusion  equation 
in  the  separator  region  is 

is£2s^^  =  ton  -  top-  (45) 

In  the  positive  electrode  region,  the  diffusive  flux  term  in  the 
right  hand  side  of  equation  (33)  gets  volume  averaged,  using  the 
boundary  conditions  (35)  and  (36)  as  follows: 


This  expression  is  the  lowest  order  approximation  which  on 
volume  averaging  will  satisfy  the  exact  result  equation  (41 ).  Higher 
order  approximations  are  not  warranted  when  uniform  reaction 
rate  approximations  are  used  for  the  reaction  rates,  as  in  equation 
(16).  This  approximation  illustrates  the  basic  approach  of  this  work: 
approximations  are  always  constructed  such  that  on  volume  aver¬ 
aging  they  respect  the  exact  results  from  volume  averaged  field 
equations.  This  approach  ensures  that  errors  vanish  at  the  domain 
(electrode  or  separator  regions)  volume  level.  The  uniform  reaction 


V.  Senthil  Kumar  /  Journal  of  Power  Sources  222  (2013)  426-441 


431 


rate  approximations,  equations  (16)  and  (17),  also  followed  this 
approach,  trivially  though. 

The  profile  approximation  (50)  is  integrated  once  with  respect 
toxas 

-  -T7  “  <*) 

where  J(t )  is  some  arbitrary  function.  Applying  the  boundary 
condition  at  separator  interface  (27)  gives 

@X=ln ■■  -q2in=-^jnin+/(t)  ^/(t)=0^D2n^=-^X.  (52) 


Note  that  equation  (52)  satisfies  the  collector  end  boundary 
condition  (25). 

Rearranging  equation  (52)  and  integrating  once  more  with 
respect  to  x  gives 


8C2 

9x 


<?2in 
lnD2  n 


c2(x,t) 


<?2ir 

kPl 


(53) 


where  J{t)  is  some  arbitrary  function.  Applying  the  interfacial 
concentration  boundary  condition  (26)  gives 

@X  =  ln:  c2in  =  -^g-|+/(t).  (54) 


Note  that  the  interfacial  concentration  appears  as  a  natural  variable. 
Hence,  it  is  also  considered  as  a  physically  relevant  internal  variable 
along  with  the  interfacial  flux.  Eliminating /(t)  between  equations 
(53)  and  (54)  gives 

ci(x,  t)  =  c2in(t)  +  (l2a  -  x2) .  (55) 

Once  the  interfacial  concentration  and  flux  are  known,  the  electro¬ 
lyte  concentration  profile  is  computed  using  equation  (55).  The  rest 
of  the  model  development  focuses  on  deriving  equations  for  the 
interfacial  concentration  and  flux.  Applying  the  volume  averaging 
definition  (7)  for  the  quadratic  term  appearing  in  equation  (55)  gives 


The  derivation  above  assumes  that  electrolyte  diffusivity  does 
not  vary  with  electrolyte  concentration,  as  per  the  cell  design 
parameters  listed  in  supplementary  material  S4.  The  derivation, 
however,  can  be  extended  to  accommodate  concentration  depen¬ 
dent  electrolyte  diffusivity.  For  example,  in  the  next  section 
concentration  dependent  ionic  conductivity  is  considered. 

For  later  reference,  from  equation  (55)  the  electrolyte  concen¬ 
tration  near  the  negative  collector  end  is  obtained  as 

c2(x  =  0,  t)  =  c2in(t)  +  .  (60) 

3.3.2.  Positive  electrode  region 

The  derivation  is  analogous  to  the  one  above,  so  only  the  salient 
steps  are  presented.  From  equation  (46),  the  profile  approximation 
is  derived  as 


This  approximation  is  integrated  with  respect  to  x  once  and  the 
boundary  condition  at  the  collector  end  (36)  is  applied  to  get 

D2^=~qf^-  (62) 


Note  that  this  equation  satisfies  the  interfacial  flux  boundary 
condition  (35). 

Integrating  once  again  with  respect  to  x  and  applying  the 
interfacial  concentration  boundary  condition  (34)  leads  to 

c!(x.t)  =  c2lp(t)-|Eg[i;-(l-^].  (63) 

Using  the  volume  average  definition  (8),  the  volume  average  of  the 
quadratic  appearing  in  the  above  equation  is 

(<‘-*)2>p4  (64> 


<*>. 


ul_il 

in  3  3  ' 


(56) 


Using  this  result  and  volume  averaging  equation  (63)  gives 
te!„«CaP-^.  (65) 


Using  this  result  and  volume  averaging  equation  (55)  gives 


(c2>n  =  c2in  + 


foin  2 
2J„D2n3 


l2 


=  c2in  + 


in<?2in 

3D2n- 


(57) 


Differentiating  this  equation  with  respect  to  time  gives 


Differentiating  this  with  respect  to  time  and  using  it  in  equation 
(47)  gives 


d(c2)n  _  dc2in  ln  dq2in 
it  dt  3D2n  if  * 


Using  the  above  expression  in  equation  (43)  gives 


,  dc2jn 
,n£2n“d T 


in£2n  dq2in 
3 D2n  dt 


— <j2in  +  (l  -t+)anln(j„). 


For  later  reference,  from  equation  (63)  the  electrolyte  concen- 

(58)  tration  near  the  positive  collector  end  is  obtained  as 

c2(x  =  L,t)  =  c2ip(t)-^^.  (67) 

(59)  3.3.3.  Separator  region 

From  equation  (44)  the  profile  approximation  is  derived  as 


Equation  (59)  is  one  equation  for  the  two  unknowns:  interfacial 
concentration  and  interfacial  flux  in  the  negative  electrode  region. 
One  more  equation  is  needed  for  these  two  unknowns  for  closure. 
These  variables  are  related  to  the  interfacial  concentration  and  flux 
in  the  positive  electrode  region,  since  an  electrolyte  balance  (48) 
prevails. 


The  electrolyte  diffusion  equation  (28)  shows  that  the  above 
expression  is  the  lowest  order  approximation  which  can  support  an 
accumulation  of  electrolyte  in  the  separator  region.  Electrolyte 
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accumulation  or  depletion  in  the  separator  is  expected  during  high 
current  transients.  This  expression  is  integrated  once  with  respect 
tox,  and  the  negative  electrode  side  flux  boundary  condition  (30)  is 
applied  to  give 

<69) 


Note  that  this  equation  satisfies  the  positive  electrode  side  flux 
boundary  condition  (32).  Integrating  once  again  with  respect  to  x 
and  applying  the  negative  electrode  side  interfacial  concentration 
boundary  condition  (29)  gives 


(x-M2 

2 


Note  that  the  positive  electrode  side  interfacial  concentration 
boundary  condition  (31)  is  not  yet  utilized.  Applying  it  on  equation 
(70)  gives  an  algebraic  equation  between  the  interfacial  concen¬ 
trations  and  fluxes 


c2in  -  C2ip 


fs(g2in  +  <J2ip) 

2D2s 


(71) 


This  equation  can  be  rearranged  as 

q  c2in  ~  c2ip  _  Q2in  +  ?2ip  {72) 


and  rationalized  as  the  relationship  between  average  concentration 
gradient  and  the  average  flux,  under  the  assumption  of  quadratic 
spatial  dependence  of  electrolyte  concentration. 

Using  the  definition  of  volume  averaging  in  the  separator  region 
(38),  the  volume  average  of  linear  and  quadratic  terms  appearing  in 
concentration  profile  (70)  are  obtained  as 


<(X-/n))s 


f|_  Is 

~  ls  2  ~  2’ 


If  s 

Is  3  3 


(73) 


Using  these  results,  volume  averaging  equation  (70)  gives 


<C2>s 


-  k<?2in  ls<hip 

2,n  3  D2s  6D2s  ' 


(74) 


Since  the  time  derivatives  of  the  average  concentrations  are  not 
independent  of  each  other,  due  to  equation  (49),  the  above  equation 
needs  not  be  differentiated  with  respect  to  time  and  used  in  the 
volume  averaged  diffusion  equation  (45)  of  the  separator.  Thus,  the 
profile  approximations  in  the  three  regions  have  yielded  three  equa¬ 
tions  (59),  (66)  and  (71)  for  the  four  unknowns:  two  interfacial 
concentrations  and  two  interfacial  fluxes.  The  fourth  equation  is  stip¬ 
ulated  among  these  four  unknowns  by  the  overall  electrolyte  balance. 
Substituting  the  average  concentrations  of  the  three  regions  (57),  (65) 
and  (74)  in  the  overall  electrolyte  balance  (48),  and  rearranging  gives 

(fn£2n  +  fs£2s)c2in+  “kin  +  £2pc2ip 

“  (30^ + 6D2SS)  92ip  =  ((n£2n + is"2s  +  /p  E2P)  C2°' 

(75) 


c2mid-c2^n  +  2,^  -  c2in  gD^  fen  8D2s  ^2ip' 


3.3.4.  Reduction  from  a  (4  x  4)  system  to  a  (2  x  2)  system 

Among  the  four  equations  for  the  four  unknowns,  two  are  ODEs 
(59)  and  (66)  while  two  are  linear  algebraic  equations  (71 )  and  (75). 
The  linear  algebraic  equations  can  be  solved  for  the  interfacial 
concentrations  in  terms  of  the  interfacial  fluxes,  and  eventually  the 
ODEs  solved  only  for  the  interfacial  fluxes.  By  this  reduction, 
instead  of  solving  4  equations  for  4  unknowns,  only  2  ODEs  are 
solved  for  the  two  interfacial  fluxes. 

Equation  (71)  is  rearranged  as 


C2  in 


Is  (<l2 in  +  <J2ip) 

C2ip  +  2D2s 


(77) 


Substituting  this  expression  for  C2m  in  equation  (75)  and  rear¬ 
ranging  gives 


c2ip  =  C20  +  ainq2in  +  aipq2ip,  (78) 

where 


l2se2s  &2n)  _ 1 _  (?9) 

,  2D2s  +  6 D2s  +  3D2n J  (/nc2n  +  lsc2s  +  /pf2p)  ’  1  J 


aip 


lnkc2n  ,  Q£2s  _  ^P£2p^  _ \ _ 

2D2s  3 D2s  3D2p J  ((nr2n  +  h^s  +  ^p£2p) 


(80) 


Differentiating  equation  (78)  with  respect  to  time  gives 

dc2ip  _  dq2in  i  dq2ip 
dt  dt  dt  ‘ 


(81) 


Similarly  differentiating  equation  (75)  with  respect  to  time,  using 
equation  (81)  in  it  for  the  time  derivative  of  c2iP  and  rearranging 
gives 


Equations  (81)  and  (82)  express  interfacial  concentration  deriva¬ 
tives  in  terms  of  interfacial  flux  derivatives.  Using  these  equations, 
interfacial  concentration  derivatives  are  eliminated  from  equations 
(59)  and  (66)  to  give 


^n£2nffin 


Un£2n  ,  ;n£2nA  dq2in  f.  Inkin']  d<72jp 

2 D2s  +  3D2n  )  dt  +  Vn  ,P  +  2 D2s  )  dt 


=  -ftan  +  0  -t+)anIn(in), 
(83) 


,  dq2in  ,  (.  ip2p\  dq2ip 

Jp^in-gp+^p^ip-^j  V 

=  <J2ip  +  (1  -  t+)ap/p(jp).  (84) 


Now,  with  four  equations  for  the  four  unknowns,  the  formulation  is 
closed. 

For  later  reference,  from  equation  (70)  the  electrolyte  concen¬ 
tration  in  the  middle  of  the  separator  is  obtained  as 


These  are  the  evolution  equations  for  the  interfacial  fluxes. 
Assuming  an  equilibrated  initial  state,  these  interfacial  fluxes  have 
zero  initial  condition.  Once  they  are  solved  for,  c2jp  is  obtained  using 
equation  (78),  and  c2jn  is  obtained  using  equation  (77).  In  the  right 
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hand  side  of  the  original  ODEs  (59)  and  (66)  only  the  interfacial 
fluxes  appear  and  not  the  interfacial  concentrations.  Hence,  it  was 
simpler  to  eliminate  interfacial  concentrations  using  the  algebraic 
equations  (71)  and  (75),  rather  than  trying  to  eliminate  the  inter¬ 
facial  fluxes. 

4.  Total  current  balance  equations 

4.1.  Field  equations 

The  total  current  at  any  point  in  the  cell  is  the  sum  of  solid  current 
(Ohmic  component  in  the  solid  phase),  migration  current  (electrolyte 
migration  due  to  electric  potential  gradient)  and  diffusion  current 
(current  carried  by  diffusive  flux  caused  by  concentration  gradient).  In 
the  separator  only  migration  and  diffusion  currents  exist  The  elec¬ 
trolyte  potential  field  is  continuous  across  the  three  regions;  hence, 
like  electrolyte  concentration  field,  electrolyte  potential  and  migration 
current  continuities  apply  at  the  electrode— separator  interfaces.  At  the 
collector  ends,  migration  currents  go  to  zero  and  solid  currents  equal 
the  external  current.  At  the  electrode- separator  interfaces,  solid 
currents  go  to  zero,  so  that  the  sum  of  migration  and  diffusion 
currents  equals  the  total  current.  The  interfacial  electrolyte  potential 
and  interfacial  migration  currents  are  identified  as  separate 
variables,  so  that  is  easy  to  impose  them  in  the  individual  regions, 
and  maintain  their  equality  across  the  interfaces. 

4.1.1.  Negative  electrode  region 


=  0.  (94) 

Having  an  electrolyte  potential  reference  within  the  separator  eases 
out  analytical  integration  in  the  reduced  order  model  development 
as  shown  below,  since  solid  current  component  is  absent  in  the 
separator  total  current  balance  (89).  The  cell  voltage  is  the  differ¬ 
ence  in  solid  potential  between  the  two  collector  ends.  Hence,  the 
choice  of  reference  potential  (both  location  and  value)  does  not 
affect  the  cell  voltage. 

4.1.3.  Positive  electrode  region 


6$i  8<2>2  ^  RT  .  9ln  c2 

+  2x2p^(l  -  4)-**  -  1, 

(95) 

with  boundary  conditions 

<f>2(/n  +  /s,t)  =  $2ip(t), 

(96) 

-«2p^  =  hip(t), 

°X  lx=i„+fs 

(97) 

H 

© 

(98) 

-  *2nl-  +  2K2n^  (1  -  t+)  - 


with  boundary  conditions 

a02|  _  n 


(85) 


4.2.  Electrolyte  potential  profiles 
For  algebraic  ease  let 


(86) 


RT 

T 


(1-4). 


(99) 


$2 On,  t)  =  $2in(t), 

(87) 

and  -K2n^  5=  hin(t)- 

(88) 

4.1.2.  Separator  region 

d&2  n  RT „  .  9ln  c2 
~^W+2^P  0-W-5T-'. 

(89) 

with  boundary  conditions 

4.2.1.  Separator  region 

The  electrolyte  concentration  profile  in  the  separator  is  known 
from  equation  (70),  using  this  profile,  equation  (89)  can  be  directly 
integrated  to  get  the  electrolyte  potential  profile  in  the  separator. 
Ionic  conductivity  is  often  a  strong  function  of  electrolyte  concen¬ 
tration.  In  the  reduced  order  framework,  however,  it  is  laborious  to 
handle  spatial  variation  of  material  properties.  Due  to  the  thinness 
of  the  electrode  or  separator  regions,  ionic  conductivity  is  assumed 
to  be  a  function  of  the  volume  average  concentration,  instead  of  the 
local  concentration: 

K2s(x,t)  =  K2[c2{x,t)\  X  £2sUg  =•  K2s(t)=K2[(C2)s(t)\  X 

(100) 


<f>2(*n,t)  =  $2in(t), 

(90) 

Note  that  the  temporal  variations  of  the  ionic  conductivity  are 

accounted  through  average  electrolyte  concentration. 

Equation 

0<2>2 1 

_,c2si#L=/n 

=  hin(t), 

(91) 

(89)  is  rearranged  to  get 

0®2  _  2@  9ln  c2  I(t) 

(101) 

<2>2(/n  +  /s,t)  • 

=  *2ip(t), 

(92) 

dx  0X  K2s(t)  ' 

,  0$2| 
and-x2s— | 

Integrating  the  above  equation  with  respect  to  x  gives 

+  ^  '2ip(t)' 

(93) 

=  2@ln  c2--j^jX+f(t), 

(102) 

One  needs 

to  assume  a  reference  potential  (either  for  electrolyte 

or  solid  potential)  somewhere  within  the  cell.  This  work  takes 

where  f[t)  is  some  arbitrary  function.  Applying  the 

potential 

a  zero  reference  potential  at  the  center  of  the  separator 

reference  (94)  gives 
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0  =  2©ln[c2mid(t)]  -  ^|y  (in  +|)  +m,  (103) 


where  the  mid-separator  concentration  is  got  from  equation  (76). 
Eliminating  J[t)  between  equations  (102)  and  (103)  gives 


$2  (x,  t)  =  2©ln 


[c2(x,t)l 

Lc2mid(t)J 


(104) 


This  is  an  exact  result,  subject  to  the  diffusion  approximations  in 
the  electrolyte  concentration  field  and  the  approximation  in 
equation  (100).  Using  the  above  expression  the  electrolyte  potential 
at  the  two  interfaces  is  identified  as 


(105) 


and *,„(,) -»2(I„  +  M)  ~^ry 

(106) 

These  interfacial  electrolyte  potentials  will  be  used  while  solving 
for  the  electrolyte  potential  in  the  electrode  regions,  to  ensure 
electrolyte  potential  continuity. 


4.2.2.  Negative  electrode  region 

As  before  the  ionic  conductivity  is  assumed  to  be  a  function  of 
the  average  electrolyte  concentration  in  this  region 

K2nW*«2[(C2>n(t)]xc^g.  (107) 


The  solid  current  term  appears  in  the  total  current  balance  for 
negative  electrode  (85).  The  expression  for  solid  current  is  taken 
from  equation  (22),  rewritten  in  a  once  integrated  form  to  facilitate 
integration  of  the  total  current  balance  equation, 


9®i  J(t)  ,, 

-'7i"DT  =  ir(,n-x) 


8  [  I(t)(ln-X)2] 

3x  [  lD  2  \ 


(108) 


Note  that  the  uniform  reaction  rate  approximation  is  entering  into 
electrolyte  potential  calculation  through  this  solid  current  equa¬ 
tion.  Using  equation  (108)  in  the  total  current  balance  (85)  gives 


8  [  f(t)  (fn-X)2]  8^2 

ax  [  "V  2  J  2n  ax 


(109) 


This  equation  can  be  rearranged  into  a  total  differential  as 
^  -  *2n*2  +  2K2n®ln  C2]  =  I.  (110) 


Ionic  conductivity,  according  to  equation  (107),  does  not  have 
spatial  dependence.  Hence,  equation  (107)  is  easily  put  in  total 
differential  form.  Integrating  once  with  respect  to  x  gives 

K2nd> 2  +  2/c2n01n  c2  =  lx  (111) 

where  J[t )  is  some  arbitrary  function.  Applying  the  interfacial 
electrolyte  potential  boundary  condition  (87)  gives 

At  x  =  /„:  -  K2n02in  +  2/c2n01n  c2in  =  I(t)ln  +f(t).  (112) 

Eliminating  f[t)  between  equations  (111)  and  (112)  gives 


02(x,t)  =  $2in(t)+20ht(|^)  +£Vx) 


1(f) 

2  K2Jn 


(fn-X)2. 


(113) 


This  is  an  exact  result,  subject  to  the  diffusion  approximations  in 
the  electrolyte  concentration  field  and  uniform  reaction  rate 
approximation  in  the  solid  current  expression. 

For  later  reference,  from  equation  (113)  the  electrolyte  potential 
near  the  negative  collector  end  is  obtained  as 


4>2(x=O,t)  =  <f>2in(t)+201n  (C?|r^) 


,  min 

2K2n 


(114) 


4.2.3.  Positive  electrode  region 

The  solution  procedure  is  similar  to  the  one  above.  The  solid 
current  equation  (23)  is  rewritten  as 

ms) 

This  expression  is  substituted  in  the  total  current  balance  (95), 
written  in  total  differential  form,  integrated  once  and  the  interfacial 
potential  boundary  condition  (96)  is  applied  to  get 


^2(x,t)  =  <P2ip(t)+201n^|g| 
|  J(t)  [X-(fn+fs)]2 

x2p/p  2 


M 

*2p 


lX-(ln+ls)} 


(116) 


For  later  reference,  from  equation  (116)  the  electrolyte  potential 
near  the  positive  collector  end  is  obtained  as 


-'-s£  <"7> 


With  equations  (104),  (113)  and  (116)  the  electrolyte  potential 
field  is  fully  solved.  The  total  current  balances  had  just  first  order 
spatial  derivatives.  Hence,  on  spatial  integration,  electrolyte 
potential  continuities  were  enough  to  fix  their  integration 
constants.  Thus,  the  migration  current  continuities,  equations  (88), 
(91),  (93)  and  (97)  were  not  utilized  yet.  These  continuities  imply 
that  the  Bruggeman  factors  for  ionic  conductivity  and  ionic  diffu- 
sivity  are  identical,  as  shown  in  the  supplementary  material  S2. 


5.  Solid  phase  mass  balance  equations 

5.1.  Field  equations 

Active  materials  in  the  negative  and  positive  electrodes  are 
idealized  as  mono-disperse  spherical  particles,  with  radii  rn  and  rp 
respectively.  Inert  lithium  species  travel  in  and  out  of  these  active 
material  spheres.  Hence,  spherical  diffusion  equations  form  the 
solid  phase  mass  balances.  The  diffusive  flux  on  the  outer  surface  of 
these  spheres  matches  the  local  electrochemical  reaction  rate.  At 
the  center  of  the  sphere  diffusive  flux  is  zero,  by  radial  symmetry. 
The  radius  of  an  active  material  sphere  is  a  measure  of  the  solid 
phase  diffusion  resistance.  The  system  of  equations  given  below  is 
identical  for  the  two  electrodes. 
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with  surface  boundary  condition  =  jk, 

symmetry  boundary  condition  =  0, 

°r  lr=0 

and  initial  condition  c1k(r,  0)  =  clk0. 


(118) 


(119) 

(120) 
(121) 


electrode  are  identical;  hence  spherical  volume  average  concen¬ 
tration  is  identical  to  electrode  volume  average  concentration. 
Hence,  under  this  approximation,  equations  (125)  and  (127)  are 
identical. 

During  discharge,  in  the  negative  electrode,  when  a  neutral  solid 
phase  Li  is  released  into  the  electrolyte  as  a  Li+  ion  and  an  electron 
is  released  into  the  outer  circuit.  This  electron  recombines  with 
some  other  Li+  in  electrolyte  within  the  positive  electrode  region 
(almost  instantaneously,  neglecting  electronic  transients),  and 
forms  a  solid  phase  Li.  Thus,  if  a  solid  phase  Li  decreases  in  the 
negative  electrode,  a  solid  phase  Li  gets  formed  at  the  positive 
electrode.  Thus,  an  overall  solid  phase  lithium  balance  holds  at  any 
time 


5.2.  Volume  averaged  field  equations 

Volume  averaging  of  a  function /(r)  within  a  sphere  is  defined  as 
follows: 

J  47tr2/(r)dr  ^ 

f  =  r=°  4  - =  ^  J  r2/(r)dr.  (122) 

2n1k  k  r  =  0 

Sphere  volume  averaging  the  left  hand  side  of  equation  (118)  shows 
that  volume  average  of  the  time  derivative  is  the  time  derivative  of 
the  volume  average 


/nrin(cin)  +  Wip<Cip>  =  infmCmo  +  WipCipo  =  Constant. 

(128) 

Differentiating  the  above  equation  with  respect  to  time  gives 

^in%^+/p^ip%^  =  0.  (129) 

This  equation  shows  that  the  rates  of  change  of  average  solid  phase 
lithium  concentration  in  the  two  electrodes  are  not  independent  of 
each  other.  For  example,  if  the  average  solid  phase  concentration  in 
the  negative  electrode  is  known,  that  in  the  positive  electrode  is 
computed  from  equation  (128)  as 

<cip>  =  (hpo+^cmo  -  (cm))-  (130) 

IpElp 


dr 


dc^ 

dt 


(123) 


The  right  hand  side  of  equation  (118)  is  sphere  volume  averaged, 
using  the  boundary  conditions  (119)  and  (120)  as 


Using  equation  (127)  in  equation  (129)  gives 


The  specific  surface  areas  of  the  active  material  spheres  are 


RHS  r2^9^)  =4  [  dr  3rlk 

r29r  V  dr  J  r3  J  r2dr\  lk  dr  J  ak 


(132) 


(124) 


Thus  the  sphere  volume  averaged  solid  state  diffusion  equation  is 


dcit.jj,  _3jk 
dt  rk  ' 


(125) 


From  the  original  initial  condition  (121 ),  the  initial  condition  for  the 
above  ODE  is  obtained  as 


Using  these  definitions  in  equation  (131)  recovers  the  electro¬ 
neutrality  condition  (13)  derived  in  Section  2.  Only  electrons 
move  in  and  out  of  a  sealed  lithium  cell.  Hence,  an  overall  lithium 
balance  is  obviously  maintained.  Due  to  electro-neutrality  at  any 
time  the  total  amounts  of  solid  phase  Li  and  electrolyte  phase  Li+ 
ions  are  individually  constant.  Hence,  there  are  separate  Li  balances 
for  the  solid  and  electrolyte  phases,  equations  (128)  and  (48) 
respectively.  As  noted  earlier,  these  volume  averaged  results  are 
exact,  even  though  they  entail  loss  of  concentration  profile  or 
gradient  information. 


cifco  =  ciko-  (326) 

The  reaction  rate  and  the  sphere  average  concentration  in  the 
above  equation  are  functions  of  x,  in  general.  Volume  averaging 
across  the  electrode  volume,  using  definitions  (7)  or  (8),  gives  the 
electrode  volume  averaged  solid  state  diffusion  equation 

d(Cifc)  _  3  (ik)  .  . 

T  rk  ■  tl2/J 

An  observation  regarding  uniform  reaction  rate  approximation, 
though  it  is  not  used  in  this  subsection:  Under  uniform  reaction  rate 
approximation,  the  concentration  fields  in  all  the  spheres  within  an 


5.3.  Profile  approximations 

The  parabolic  and  quartic  profile  approximations  for  the  solid 
state  diffusion  in  a  lithium  cell  are  reported  in  Ref.  [10],  wherein  it  is 
shown  that  high  current  scenarios  warrant  a  quartic  approxima¬ 
tion.  For  the  sake  of  completeness,  consistent  style  (volume  aver¬ 
aging  based)  and  notation,  the  quartic  approximation  is  reworked 
below.  The  parabolic  approximation  derivation  is  given  in  the 
supplementary  material  S3. 

For  the  lowest  order  volume  average  respecting  parabolic 
approximation,  the  profile  approximation  was  derived  from  Equa¬ 
tion  (124)  as 
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i  a 

r*dr 


'  dr  ) 


rk  ' 


(142) 


The  right  hand  side  of  the  above  equation  is  independent  of  the 
radial  coordinate;  it  is  dependent  only  on  time  through  the  reaction 
rate.  For  a  higher  order  approximation,  the  gradient  of  diffusive  flux 
is  assumed  to  be  parabolic  in  the  radial  coordinate 

PS  0*1?)  +  t'34) 


This  is  a  second  equation  for  the  two  unknowns  G  and  H.  Equations 
(136)  and  (142)  are  solved  simultaneously  to  get  explicit  expres¬ 
sions  for  G  and  H  as 


-clfcr;  H  =  ——jk - 


Integrating  equation  (139)  once  more  with  respect  to  r  gives 


From  the  sphere  volume  average  definition  (122),  the  volume  cifc  =  cfi  P  + 
average  of  r2  is  obtained  as  6Dlk  20Dlkrj* 


(144) 


(135) 


Sphere  volume  averaging  equation  (134)  using  the  above  result 
gives 


1  f 
r^dr 


k  dr  ) 


% 

rk  ' 


(136) 


This  gives  one  equation  for  the  two  unknowns  G  and  H.  Rearranging 
equation  (134)  and  integrating  once  with  respect  to  r  gives 


0 

dr 


(r2Dlk 


dclk 

dr 


s) 


=  Gr2  +  H^  => 

-4+h^+«(>' 


(137) 


where  f{t)  is  some  arbitrary  function.  Applying  the  surface 
concentration  condition  to  the  above  equation  gives 


mhs'k+m- 

i  k  20  Dur£ 


(145) 


Eliminating  J[t)  between  equations  (144)  and  (145)  gives  the 
quartic  concentration  profile 

-2®b <146) 

For  sphere  volume  averaging  the  right  hand  side,  equation  (146) 
and  the  following  result  are  needed: 


lrl  =  3  4 

rl  7  7  k' 


(147) 


where  J(t )  is  some  arbitrary  function.  Using  the  symmetry 
boundary  condition  (120)  it  is  easily  seen  that/(t)  is  zero.  Hence, 


r2D-[k 


dcik  rr3 
G  3 


H^k- 


(138) 


Sphere  volume  averaging  equation  (146)  and  rearranging  gives 

c**s'*+llS^  +  3sfeH-  (,48> 

Substituting  the  expression  (143)  for  G  and  H  in  the  above  equation 
and  simplifying  gives 


It  is  easily  seen  that  applying  the  surface  flux  boundary  condition 
(119)  reduces  the  above  equation  to  (136).  Hence,  the  surface  flux 
boundary  condition  does  not  generate  another  equation  for  the  two 
unknowns.  Equation  (138)  is  rearranged  as 


0^ 

dr 


ft 

5Dikrfc 


The  average  radial  gradient 


(139) 


(140) 


Csk  =  ft k  -  +  S  c ifer-  (149) 

Under  uniform  reaction  rate  approximation,  jk(x,  t)~(jK)(t)  is 
known.  The  spherical  average  concentrations  are  obtained  by 
solving  the  sphere  volume  averaged  solid  diffusion  equation  (125). 
The  rest  of  the  subsection  derives  an  evolution  equation  for  the 
spherical  average  radial  gradient  c1)ir.  Once  it  is  identified,  equation 
(149)  is  used  to  feed  surface  concentrations  to  the  Butler-Volmer 
kinetics,  as  shown  in  the  next  section. 

Differentiating  the  solid  state  diffusion  equation  (118)  with 
respect  to  r  and  using  the  equality  of  mixed  derivatives  gives 


is  introduced  as  a  physically  relevant  variable;  an  evolution  equa¬ 
tion  for  it  will  be  derived  by  the  end  of  this  subsection.  For  sphere 
volume  averaging  the  linear  and  cubic  terms  in  the  right  hand  side 
of  (139)  the  following  results  are  needed: 


Using  equations  (140)  and  (141 ),  sphere  volume  averaging  equation 
(139)  and  simplifying  gives 


im  -tm  -'so, 

Using  the  profile  approximation  (134)  for  the  radial  diffusive  flux  in 
the  right  hand  side  gives 


Sphere  volume  averaging  the  above  equation,  using  equation  (141) 
for  the  averaging  of  the  linear  term  on  the  right  hand  side,  and 
simplifying  gives 
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dclkr  _  3  H 

*W’'T  2Tk 


(152) 


Substituting  for  H  from  equation  (143)  and  simplifying  gives  the 
evolution  equations  for  the  average  radial  gradients  as 


dcikr 

IT 


30  Dlk  45 

,3  c' kr  =  ~^2Jk- 


(153) 


These  ODEs  for  the  average  radial  gradient  are  solved  with  the 
initial  conditions 


clfcr(r,Q)  =  Q,  (154) 

when  the  initial  states  are  in  equilibrium,  with  no  concentration 
gradients. 


potential  field.  In  the  solid  phase,  using  the  reaction  rate  the  solid 
phase  concentration  field  is  solved.  All  the  above  fields  provide 
input  for  the  Butler-Volmer  kinetics  to  solve  for  the  solid  potential 
field.  This  computational  methodology  is  used  for  both  the 
electrodes.  Eventually  the  cell  voltage  is  computed  as  the  solid 
potential  difference  between  the  two  collector  ends. 

For  the  lowest  order  reduced  order  model  derived  here,  the 
reaction  rates  in  the  electrodes  are  assumed  to  be  the  electrode 
volume  average  reaction  rates.  The  salient  results  of  this  model  are 
presented  below,  in  a  sequential  algorithmic  order.  This  reduced 
order  model  predicts  not  just  the  cell  voltage,  but  also  the  profiles 
of  all  the  internal  variables,  subject  to  the  uniform  reaction  rate  and 
profile  based  approximations.  The  cell  voltage  and  predictions  from 
the  reduced  order  model  are  compared  with  the  detailed  electro¬ 
chemical  model  results  in  the  next  section. 


6.  Butler-Volmer  kinetics 


7.1.  Reaction  rates  and  solid  current  profiles  in  the  electrodes 


Butler-Volmer  kinetics  for  the  electrochemical  reactions 
occurring  on  the  surface  of  the  active  material  spheres  in  the 
electrodes  are 

Jk  =  2/fcoSinh  (#i  -  $2  -  14)] ;  k=  n,  p,  (155) 

where 

jko  =  MCskmax  -  CSk)°  5  <%k  <%  5  ’  (156) 

and  Uk  are  the  electrode  open  circuit  potentials. 

Inverting  equation  (155)  and  rearranging  gives 

^  =Uk  +  <P2+  ^sinh-1  (J^j .  (157) 

This  equation  is  used  to  find  the  solid  potential  at  any  point  within 
the  electrodes,  with  the  following  inputs  from  the  earlier  sections: 

(1)  The  reaction  rate  4  is  obtained  using  the  uniform  reaction  rate 
approximation  equations  (16)  or  (17). 

(2)  The  electrolyte  concentration  C2,  appearing  in  jko,  is  obtained 
from  equation  (55)  or  (63). 

(3)  The  electrolyte  potential  $ 2  is  obtained  from  equation  (113)  or 
(116). 

(4)  The  solid  phase  concentration  csk,  appearing  in  jk0,  is  obtained 
from  equation  (148)  for  the  quartic  approximation. 

The  cell  voltage  is  the  solid  potential  difference  between  the 
positive  and  negative  collector  ends  i.e. 

V  =  $,(x  =  I)-$i(x  =  0).  (158) 

Since  the  reduced  order  model  derivation  is  spread  out  across 
various  sections,  and  the  computational  order  is  different  from  the 
derivation  order,  the  model  results  are  collated  in  the  next  section 
as  a  ready  to  use  algorithm. 

7.  Reduced  order  model  with  uniform  reaction  rate 
approximation  -  algorithm 

In  the  reduced  order  model  framework,  with  a  consistent  reac¬ 
tion  rate  approximation,  the  electrolyte  and  the  solid  phase  fields 
are  solved  in  a  particular  order  as  shown  in  Fig.  2.  In  the  electrolyte 
phase,  using  the  reaction  rate  the  electrolyte  concentration  field  is 
solved  first,  which  then  enables  the  solution  of  the  electrolyte 


External  current  fixes  the  average  reaction  rate  in  the  electrodes, 
at  any  time.  Under  the  uniform  reaction  rate  approximation,  all 
points  within  an  electrode  experience  the  same  rate.  The  uniform 
reaction  rates  in  the  negative  and  positive  electrodes  are  given  by 
equations  (16)  and  (17),  respectively.  Uniform  reaction  rate 
approximation  implies  that  solid  current  profiles  are  linear.  The 
linear  current  profiles  in  the  negative  and  positive  electrodes  are 
given  by  equations  (22)  and  (23),  respectively.  These  internal 
profiles  are  not  required  for  cell  voltage  computation. 

7.2.  Electrolyte  concentration  field  in  the  electrode  and  separator 
regions 

The  electrolyte  concentration  field  is  expressed  in  terms  of  the 
electrode— separator  interfacial  concentrations  and  fluxes.  The 
interfacial  fluxes  are  obtained  by  solving  the  two  coupled  ODEs, 
equations  (83)  and  (84),  with  zero  initial  condition  for  the  inter¬ 
facial  fluxes  for  an  equilibrated  initial  state.  The  parameters  am  and 
aip  appearing  in  these  two  ODEs  are  obtained  from  equations  (79) 
and  (80),  respectively.  The  interfacial  concentrations  are 
computed  from  the  interfacial  fluxes  using  equations  (77)  and  (78). 
For  computing  the  cell  voltage  the  electrolyte  concentrations  at  the 
collector  ends  and  the  middle  of  the  separator  are  needed,  they  are 
obtained  from  equations  (60),  (67)  and  (76).  The  electrolyte 


Kinetics 


Solid  potential 
field 


Fig.  2.  Reduced  order  model  computational  methodology  for  an  electrode. 
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concentration  profiles  in  the  three  regions  are  computed  using 
equations  (55),  (63)  and  (70).  These  internal  profiles  are  not 
required  for  cell  voltage  computation. 

7.3.  Electrolyte  potential  field  in  the  electrode  and  separator  regions 

The  electrolyte  potential  field  is  expressed  in  terms  of  the 
electrode-separator  interfacial  potentials.  In  this  work  zero 
reference  potential  is  set  in  the  middle  of  the  separator.  These 
interfacial  electrolyte  potentials  are  computed  from  equations 
(105)  and  (106).  Note  that  these  equations  use  the  interfacial  and 
mid-separator  electrolyte  concentrations  computed  in  the  earlier 
subsection.  The  electrolyte  potential  at  the  collector  ends  are 
computed  using  these  interfacial  electrolyte  potentials  using  equa¬ 
tions  (114)  and  (117).  The  electrolyte  potential  profiles  in  the  three 
regions  are  computed  using  equations  (104),  (113)  and  (116).  These 
internal  profiles  are  not  required  for  cell  voltage  computation. 

7.4.  Solid  phase  concentration  field 

Quartic  profile  approximation  expresses  the  concentration 
profile  within  an  active  material  sphere  in  terms  of  the  average 
concentration  and  the  average  radial  gradient.  Under  the  uniform 
reaction  rate  approximation,  all  the  spheres  in  an  electrode  have 
identical  concentration  field  at  all  times  (since  initial  concentration 
and  reaction  rates  are  identical).  For  the  negative  electrode  the 
following  two  coupled  ODEs,  equations  (125)  and  (153)  are  solved 
with  initial  conditions  (126)  and  (154).  The  surface  concentration  of 
the  negative  electrode  spheres  is  given  by  equation  (149). 

Instead  of  solving  an  ODE  for  the  average  solid  phase  Li 
concentration  in  the  positive  electrode,  it  is  computed  from  the 
solid  phase  lithium  balance  using  equation  (130).  The  electrode 
volume  average  concentration  is  identical  to  the  spherical  volume 
average  concentration  under  uniform  reaction  rate  approximation. 
For  the  average  radial  concentration  gradient  in  the  positive  elec¬ 
trode  sphere,  the  ODE  (153)  is  solved  with  the  initial  condition 
(154).  The  surface  concentration  of  the  positive  electrode  spheres  is 
given  by  equation  (149). 

The  concentration  profiles  inside  the  active  material  spheres  are 
computed  in  terms  of  surface  concentration,  average  radial 
gradient,  and  reaction  rate  using  equations  (143)  and  (146).  These 
internal  profiles  are  not  required  for  cell  voltage  computation. 

7.5.  Butler—  Volmer  kinetics 

In  this  reduced  order  framework,  as  shown  in  Fig.  2,  reaction  rates, 
electrolyte  concentrations,  electrolyte  potentials  and  solid  phase 
concentrations  feed  into  Butler-Volmer  kinetics,  which  is  then 
analytically  inverted  to  compute  the  solid  potential.  The  reaction 
rate  pre-factor  in  the  negative  electrode  is  computed  using  equa¬ 
tion  (156).  The  electrolyte  concentration  and  solid  phase  surface 
concentrations  therein  are  obtained  from  the  earlier  subsections.  The 
solid  potential  in  the  two  electrodes  are  computed  using  equation 
(157).  The  electrolyte  potential  and  the  reaction  rate  used  therein  are 
obtained  from  the  earlier  subsections.  For  cell  voltage  computation, 
equation  (157)  is  evaluated  at  the  two  collector  ends:  x  =  0  and  x  =  L. 
Eventually,  the  cell  voltage  is  computed  as  the  difference  in  solid 
potential  at  the  end  of  the  two  collectors  as  in  equation  (158). 

7.6.  Model  reduction  summary 

The  detailed  electrochemical  isothermal  model  solves  10 
coupled  partial  differential  equations:  two  solid  phase  current 
balance  equations  (1)  and  (4),  three  electrolyte  diffusion  equations 
(24),  (28)  and  (33),  three  total  current  balance  equations  (85),  (89) 


and  (95),  and  two  solid  phase  diffusion  equations  —  equation  (118) 
applied  for  negative  and  positive  electrode  spheres  individually. 
The  isothermal  reduced  order  model  reported  herein  solves  5  linear 
differential  equations:  two  for  the  interfacial  fluxes  equations  (83) 
and  (84),  one  for  the  average  concentration  in  the  negative  elec¬ 
trode  sphere  equation  (125),  and  two  for  the  average  concentration 
gradients  in  the  negative  and  positive  electrode  spheres  equation 
(153).  Thus  instead  of  solving  10  coupled  PDEs,  this  ROM  solves 
only  5  linear  ODEs.  This  complexity  reduction  was  achieved  using 
uniform  reaction  rate  and  profile  approximations. 

This  reduced  order  model  is  structured  such  that  nonlinearities 
are  confined  to  the  final  algebraic  evaluations  (Computation  of 
electrolyte  potential  by  analytical  integration,  e.g.  Equation  (113) 
and  inversion  of  Butler-Volmer  kinetics,  e.g.  Equation  (157))  and 
the  five  ODEs  being  solved  are  linear.  The  ideal  model  reduction 
that  can  be  achieved  is  abstracted  in  Fig.  3.  A  dynamical  system  of 
engineering  interest  is  often  described  by  nonlinear  PDEs.  Through 
suitable  approximations  one  can  reduce  the  dynamical  information 
into  a  system  of  linear  ODEs  and  relegate  the  nonlinearities  to 
algebraic  expressions,  as  in  this  work.  The  nonlinearity  of  the 
problem,  once  relegated  to  algebraic  expressions,  is  not  bother¬ 
some  because  linear  coupled  ODEs  can  be  solved  without  iteration. 
Such  a  time  marching  algorithm  with  minimal  computational 
complexity  is  needed  for  on-board  control  applications. 

7.7.  Voltage  transients  during  rest 

Uniform  reaction  rate  approximation  equates  the  local  reaction 
rate  everywhere  within  an  electrode  to  the  volume  average  reaction 
rate,  and  the  volume  average  reaction  rate  is  proportional  to  the 
external  current  at  any  time.  Hence,  if  external  current  is  cut  off, 
under  this  approximation,  reaction  rate  goes  to  zero  instantaneously. 
One  of  the  causes  for  the  voltage  transients  during  rest  is  the  relax¬ 
ation  of  concentration  gradients  within  active  material  spheres.  This 
phenomenon  is  captured  by  the  quartic  approximation  for  solid 
diffusion  presented  above.  Since  the  local  reaction  rate  goes  to  zero 
instantaneously  under  uniform  reaction  rate  approximation,  the 
average  concentration  within  the  active  material  spheres  remains 
frozen  according  to  the  equation  (125).  The  surface  concentration  in 
an  active  material  sphere  given  by  equation  (149),  however,  can 
change  as  the  average  radial  gradient  decreases  during  rest.  Once  the 
surface  concentration  feeding  into  Butler-Volmer  kinetics  changes, 
solid  potential  changes,  leading  to  voltage  transients.  The  right 
hand  side  of  equation  (149)  shows  the  two  pathways  by  which 
surface  concentration  relaxes  during  rest:  one  by  reaction  rate  and 
another  by  radial  concentration  gradient.  Most  of  the  voltage  relax¬ 
ation  during  rest  is  expected  to  arise  from  the  relaxation  of 
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Fig.  3.  Ideal  model  reduction  methodology. 


V.  Senthil  Kumar  /  Journal  of  Power  Sources  222  (2013)  426-441 


439 


concentration  gradients;  this  pathway  is  captured  by  the  current 
model.  The  spatially  varying  reaction  rate  can  also  contribute  to 
voltage  relaxation  during  rest;  this  pathway  is  not  captured  by  the 
current  model. 

8.  Comparison  of  reduced  order  model  predictions  with  the 
detailed  electro-chemical  model  results 

In  this  section  cell  voltage  predictions  from  the  reduced  order 
model  with  uniform  reaction  rate  approximation  are  compared 
with  those  from  the  detailed  electrochemical  model  (default  model 
in  Comsol  Multiphysics®)  for  a  typical  commercial  cell  for  constant 
current  and  pulse  protocols.  The  cell  design  parameters  used  in  the 
Comsol®  model  are  presented  in  the  supplementary  material  S4. 

8.1.  Constant  current  protocols 

The  cell  voltage  predictions  for  1 C  constant  current  discharge  are 
shown  in  Fig.  4,  the  absolute  and  percentage  errors  in  the  ROM 
prediction  in  Figs.  5  and  6,  respectively.  At  1C  constant  current 
discharge  ROM  gives  a  maximum  error  of  about  9  mV  and 
a  maximum  percentage  error  of  about  0.27%  as  shown  by  these 
figures.  The  corresponding  figures  for  7C  constant  current  discharge 
are  Figs.  7—9,  which  show  that  the  current  ROM  gives  a  maximum 
error  of  about  80  mV  and  a  maximum  percentage  error  of  about  3.2%. 
It  is  interesting  to  note  that  with  uniform  reaction  rate  approxima¬ 
tion  and  the  leading  order  quadratic  approximation  for  the  elec¬ 
trolyte  concentrations  in  the  three  regions,  the  cell  voltage 
prediction  has  less  than  3.2%  error,  even  at  high  current  discharge. 

8.2.  Pulse  protocols 

A  typical  low  current  pulse  protocol  is  shown  in  Fig.  10.  It  has 
2000  s  of  1C  discharge,  followed  by  300  s  of  rest,  2000  s  of  1C 
charge,  and  rest  till  8000  s.  The  cell  voltage  predictions  are 
compared  in  Fig.  11,  and  the  percentage  error  in  ROM  prediction  is 
shown  in  Fig.  12.  The  latter  figure  shows  that  the  maximum 
percentage  error  is  about  0.2%.  The  maximum  absolute  error  (figure 
not  shown)  is  about  10  mV,  occurring  during  constant  current 
periods. 

A  typical  high  current  pulse  protocol  is  shown  in  Fig.  13.  It  has 
10  s  of  7C  discharge  followed  by  40  s  of  rest,  10  s  of  7C  charge,  and 
60  s  of  rest.  The  cell  voltage  predictions  are  compared  in  Fig.  14,  and 
the  percentage  error  in  ROM  prediction  is  shown  in  Fig.  15.  The 
maximum  percentage  error  is  about  2%  as  shown  in  the  latter 
figure.  The  maximum  absolute  error  (figure  not  shown)  is  about 
80  mV,  occurring  during  constant  current  periods. 


Fig.  5.  Absolute  error  in  the  cell  voltage  predicted  by  the  reduced  order  model  at  1C 
constant  current  discharge. 

It  is  interesting  to  note  that  the  errors  in  the  cell  voltage 
predictions  in  both  these  pulse  protocols  during  the  rest  periods  are 
negligible  (in  Fig.  12  period  2000-2300  s,  in  Fig.  15  periods  10-50  s 
and  60-120  s).  These  results  show  that  the  voltage  transients  at  the 
beginning  of  the  rest  periods  are  well  captured  by  ROM  at  the 
current  levels  of  approximation. 

9.  Summary 

The  detailed  electrochemical  model  for  a  single  cell  comprises  of 
10  coupled  partial  differential  equations  (PDEs),  for  an  isothermal 
operation.  It  is  computationally  intensive  to  use  this  model  for 
battery  packs,  battery  management  system  algorithms,  calendar  or 
cycle  life  predictions,  detailed  parameter  estimation,  inclusion  of 
complex  degradation  mechanisms  etc.  This  work  describes 
a  reduced  order  model  (ROM)  framework,  which  reduces  the 
detailed  model  PDEs,  to  a  manageable  set  of  ordinary  differential 
equations  (ODEs),  which  can  be  used  for  the  above  applications. 

Volume  averaging  PDEs  yield  ODEs.  Hence  volume  averaging  is 
the  basic  method  used  for  model  reduction.  The  gradient  or  profile 
information  lost  on  volume  averaging  is  recovered  through  profile 
based  approximations.  The  electrolyte  concentration  profiles  are 
expressed  in  terms  of  electrode— separator  interfacial 
concentrations  and  fluxes.  From  the  detailed  model,  equations  are 
derived  for  these  physically  relevant  internal  variables. 

In  the  ROM  framework,  with  a  consistent  reaction  rate  approxi¬ 
mation,  the  electrolyte  and  the  solid  phase  fields  are  solved  as 
follows;  In  the  electrolyte  phase,  using  the  reaction  rate  the  elec¬ 
trolyte  concentration  field  is  solved  first,  which  then  enables  the 
solution  of  the  electrolyte  potential  field.  In  the  solid  phase,  using  the 
reaction  rate  the  solid  phase  concentration  field  is  solved.  All  the 
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Fig.  6.  Percentage  error  in  the  cell  voltage  predicted  by  the  reduced  order  i 
constant  current  discharge. 


:  1C 


440 


V.  Senthil  Kumar  /  Journal  of  Power  Sources  222  (2013)  426—441 


Fig.  7.  Cell  voltage  predictions  from  the  detailed  electrochemical  model  (Comsol®)  and  FiS-  «•  Cell  voltage  predictions  from  the  detailed  electrochemical  model  (Comsol®) 
the  reduced  order  model  (ROM)  at  7C  constant  current  discharge.  and  the  reduced  order  model  (ROM)  for  the  low  current  pulse  protocol  in  Fig.  10. 


Fig.  8.  Absolute  error  in  the  cell  voltage  predicted  by  the  reduced  order  model  at  7C 
constant  current  discharge. 
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Fig.  9.  Percentage  error  in  the  cell  voltage  predicted  by  the  reduced  order  model  at  7C 
constant  current  discharge. 


Fig.  12.  Percentage  error  in  the  cell  voltage  predicted  by  the  reduced  order  model  for 
the  low  current  pulse  protocol  in  Fig.  10. 


Fig.  13.  A  typical  high  current  pulse  protocol  used  to  compare  the  reduced  order 
model  and  the  detailed  electrochemical  model. 


-20 


0 


2000  4000  6000  8000 

t(s) 


t(s) 


Fig.  10.  A  low  current  pulse  protocol  used  to  compare  the  reduced  order  model  and  Fig.  14.  Cell  voltage  predictions  from  the  detailed  electrochemical  model  (Comsol®) 
the  detailed  electrochemical  model.  and  the  reduced  order  model  (ROM)  for  the  high  current  pulse  protocol  in  Fig.  13. 
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Fig.  15.  Percentage  error  in  the  cell  voltage  predicted  by  the  reduced  order  model  for 
the  high  current  pulse  protocol  in  Fig.  13. 

above  fields  provide  input  for  the  Butler-Volmer  kinetics  to  solve  for 
the  solid  potential  field.  This  computational  methodology  is  used  for 
both  the  electrodes.  Eventually  the  cell  voltage  is  computed  as  the 
solid  potential  difference  between  the  two  collector  ends. 

For  the  lowest  order  ROM,  a  uniform  reaction  rate,  quadratic 
electrolyte  concentrations  and  quartic  solid  phase  concentrations 
are  assumed.  In  this  ROM  only  5  linear  ODEs  are  solved,  and  the  rest 
are  nonlinear  algebraic  expressions.  When  compared  with  the 
detailed  electrochemical  model,  this  ROM  gives  a  maximum  error 
of  about  0.27%  during  1C  discharge  and  about  3.2%  during  7C 
discharge,  for  a  typical  commercial  cell. 
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Notation 


Constants  and  variables 
A:  electrode  cross  sectional  area,  m2 

ak  =  3r-ik/rk:  specific  surface  area  of  active  materials  in  n  and  p  electrodes,  m-1 
brug:  Bruggeman  factor,  used  to  calculate  effective  porous  media  properties, 
dimensionless 

Cj,  a.  solid  phase  concentration  in  active  material  spheres  and  electrolyte  concen- 

csfc  Cj*  max:  surface  solid  phase  concentrations,  and  their  maximum  value,  mol  m-3 
C2ik,  Q2ik:  interfacial  electrolyte  concentrations  and  fluxes,  mol  m-3,  mol  m-2  s_1 
D2:  electrolyte  diffusivity  (material  property),  m2  s_1 

D2 k  =  D2  x  £2™S,;  effect've  electrolyte  diffusivity  in  the  porous  k  =  n,s,p  domains 
Du,:  solid  diffusivity  of  k  =  n,p  electrode  material  spheres,  m2  s_1 
F:  Faraday's  constant,  96.487C  mol-1 
1:  external  current  density,  A  m-2 

ji ,:  local  surface  reaction  rate  in  k  =  n,  p  electrodes,  mol  m-2  s_1 

ik°  mol  m-^“  -  csJa5c?jfc5cf5.’  over-potential  independent  rate  pre-factor, 

kfc;  surface  reaction  rate  constant  in  k  =  n,  p  electrodes,  (mol  m-2  s_1)(mol  m_3)_l  s 
In,  Is,  Ip:  thickness  of  negative  electrode,  separator  and  positive  electrode,  m. 

L  =  ln  +  (s  +  lp  is  total  cell  thickness 
R:  gas  constant,  8.314  J  mol-1  K_1 

r„,  rp:  Radii  of  active  material  spheres  in  the  negative  and  positive  electrodes,  m 

SOCk  =  Csk/Cskm ax-'  state  of  charge  of  k  =  n,  p  electrodes,  dimensionless 

T:  temperature,  298  K  in  this  work 

t+:  electrolyte  transference  number,  dimensionless 

Uk:  open  circuit  voltage  of  k  =  n,  p  electrodes,  V 

£i  k,  £2k,  £fk:  volume  fractions  of  active  material,  electrolyte  and  filler  in  electrode  or 
separator  domains,  dimensionless.  eu*  +  C2k  +  £ot  =  1.  Separator  has  no  active 
material.  Hence  £is  =  0 

k2:  electrolyte  conductivity  (material  property),  S  m-1 

k-2 k  =  k2x  £2™s:  effective  electrolyte  conductivity  in  porous  k  =  n,s,p  domains 
a\c.  electrical  conductivity  of  k  =  n,p  electrode  materials  (material  property),  S  m_1 
ff] k  =  <r*  x  cb™g:  effective  electrical  conductivity  in  the  porous  k  =  n,p  domains, 
S  m_1 


Supplementary  material  associated  with  this  article  can  be 
found,  in  the  online  version,  at  http://dx.doi.org/10.1016/jjpowsour. 
2012.09.013. 
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$i,  #2-'  solid  and  liquid  potential,  V 
Operators  and  subscripts 

(■):  volume  average,  defined  in  the  three  domains,  n,s,p 

•:  volume  average  within  an  active  material  sphere  in  n  or  p  electrode 

J,  2,  f:  solid  active  material,  liquid  electrolyte  phases  and  solid  fillers  (including 
additives) 
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in:  separator  interface  of  the  negative  electrode 
ip:  separator  interface  of  the  positive  electrode 
k  =  n,s,p:  negative  electrode,  separator,  and  positive  electrode 
r:  radial  gradient 


